Regression (Period-based) CLEAN

Intro

Here we provide exploration of LMER models fitted to the 2020 U.S. Congressional election Twitter Dataset (https://doi.org/10.6084/m9.figshare.25523734).

This document was typeset assembled using Quarto (https://quarto.org).

Initialize Workspace

Clear and Init Workspace

Imports

library(ggplot2)
library(GGally)
library(plyr)
library(tidyverse)

Set Directories

Data Load and Inspection

Data Load

ifn0 <- "i0001-politicians-aggregated.rds"
suppressWarnings(rm(list = ls(pattern = "^df")))
df0 <- readr::read_rds(file.path(ifd0, ifn0)) %>%
  dplyr::select(Name, Party, Outcome, Period, Time, Days, Agency) %>%
  dplyr::mutate(
    Period = forcats::fct_recode(
      Period,
      PC = "MI",
      PB = "BE",
      PA = "AR")) %>%
  identity()
df0 %>% dim() %>% cat0()

169997 7 

for (name in names(df0)) {cat0(paste0("    ", name, ", "))}

    Name,  
    Party,  
    Outcome,  
    Period,  
    Time,  
    Days,  
    Agency,  

Skim

set.w(222)
library(skimr)
options(scipen=999)
skimr::skim_format(numeric=list(digits=6))

Warning: Formatting options are now in print() or as function arguments in skim_with().

## digits=4, scientific = FALSE
df0 %>% skimr::skim() %>% dplyr::select(-c(n_missing)) %>% print(numeric=list(digits=2))

── Data Summary ────────────────────────
                           Values    
Name                       Piped data
Number of rows             169997    
Number of columns          7         
_______________________              
Column type frequency:               
  factor                   4         
  numeric                  3         
________________________             
Group variables            None      

── Variable type: factor ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate ordered n_unique top_counts                            
1 Name                      1 FALSE        870 Bar: 361, Ela: 361, Llo: 361, Pra: 361
2 Party                     1 FALSE          2 Dem: 99081, Rep: 70916                
3 Outcome                   1 FALSE          2 win: 115031, los: 54966               
4 Period                    1 FALSE          3 PB: 98660, PA: 47974, PC: 23363       

── Variable type: numeric ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
  skim_variable complete_rate     mean      sd      p0     p25     p50    p75   p100 hist 
1 Time                      1  -0.0680   0.572   -1     -0.55   -0.122  0.428   1    ▇▇▇▆▆
2 Days                      1 -12.2    103.    -180    -99     -22     77     180    ▇▇▇▆▆
3 Agency                    1   0.500    0.271   -1.69   0.343   0.500  0.658   2.37 ▁▁▇▂▁

Check factors order

df2 <- df0
df3 <- df0
summary(df0$Period)

   PC    PB    PA 
23363 98660 47974 

df2$Period <- factor(df2$Period, levels = c("PB", "PC", "PA"))
summary(df2$Period)

   PB    PC    PA 
98660 23363 47974 

df3$Period <- relevel(df0$Period, ref = "PB")
summary(df3$Period)

   PB    PC    PA 
98660 23363 47974 

Prepare for Regression

File names and output directories

## ./data/j1002-model-by-period/i0001_models/fit0x0.extension

Stuff for plotting

Marginalization

Checkups

Save data frame for reference

file <- file.path(dir8, "data-df0.rds")
df0 %>% readr::write_rds(file)
cat0(file)

./data/j1002-model-by-period/i0001_models/data-df0.rds 

Model 01: Null

Fit

model <- "fit01aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (1 | Name) + 1,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit01aPd: [df0] Agency ~ (1 | Name) + 1

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (1 | Name) + 1
   Data: df0
Control: control

REML criterion at convergence: 26631.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-7.9514 -0.5639 -0.0023  0.5683  7.3986 

Random effects:
 Groups   Name        Variance Std.Dev.
 Name     (Intercept) 0.006209 0.0788  
 Residual             0.067519 0.2598  
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value            Pr(>|t|)    
(Intercept)   0.498618   0.002807 826.426886   177.7 <0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit01aPd: [df0] Agency ~ (1 | Name) + 1

# R2 for Mixed Models

  Conditional R2: 0.084
     Marginal R2: 0.000

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit01aPd: [df0] Agency ~ (1 | Name) + 1

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.084
  Unadjusted ICC: 0.084

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit01aPd: [df0] Agency ~ (1 | Name) + 1

# ICC by Group

Group |   ICC
-------------
Name  | 0.084

Random Effects

extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 02: Time

Fit

model <- "fit02aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit02aPd: [df0] Agency ~ (Time | Name) + Time

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time
   Data: df0
Control: control

REML criterion at convergence: 24941.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.0353 -0.5639 -0.0047  0.5664  7.3362 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006283 0.07927      
          Time        0.004088 0.06394  0.18
 Residual             0.066384 0.25765      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
              Estimate Std. Error         df t value             Pr(>|t|)    
(Intercept)   0.498029   0.002848 816.443563 174.884 < 0.0000000000000002 ***
Time         -0.010230   0.002634 709.116999  -3.884             0.000112 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr)
Time 0.187 

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit02aPd: [df0] Agency ~ (Time | Name) + Time

# R2 for Mixed Models

  Conditional R2: 0.102
     Marginal R2: 0.000

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit02aPd: [df0] Agency ~ (Time | Name) + Time

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.102
  Unadjusted ICC: 0.102

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit02aPd: [df0] Agency ~ (Time | Name) + Time

Model contains random slopes. Cannot compute accurate ICCs by group factors.

# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

## pp(model);
tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)
tbl0

Characteristic	Beta	95% CI1	p-value
Time	-0.01	-0.02, -0.01	<0.001
Name.sd__(Intercept)	0.08
Name.cor__(Intercept).Time	0.18
Name.sd__Time	0.06
Residual.sd__Observation	0.26
1 CI = Confidence Interval

Model Parameters

pp(model);

fit02aPd: [df0] Agency ~ (Time | Name) + Time

parameters::model_parameters(get(model))

# Fixed Effects

Parameter   | Coefficient |       SE |         95% CI | t(169991) |      p
--------------------------------------------------------------------------
(Intercept) |        0.50 | 2.85e-03 | [ 0.49,  0.50] |    174.88 | < .001
Time        |       -0.01 | 2.63e-03 | [-0.02, -0.01] |     -3.88 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est") 
gg88 <- gg88 + ggplot2::ylim(-0.1,0.1) + line2 + stuff0

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time

terms <- c("Time")
extra <- "0001"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + time0 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

pp(model);

fit02aPd: [df0] Agency ~ (Time | Name) + Time

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.51 | 0.50, 0.51
-0.50 |      0.50 | 0.50, 0.51
 0.00 |      0.50 | 0.49, 0.50
 0.50 |      0.50 | 0.49, 0.50
 1.00 |      0.49 | 0.48, 0.50

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Slope     |       95% CI |     p
--------------------------------
-7.57e-03 | -0.01,  0.00 | 0.004

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.15,)

gg88 <- gg88 + line0 + line1 + time0 + stuff0
## gg88 <- gg88 + coord_cartesian(ylim = c(0.45, 0.55))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

Random Effects

extra <- "9001"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="re")

file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-random-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=4, height=88, limitsize = FALSE)
gg88

Model 03: Time x Period

Fit

model <- "fit03aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Period,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Period
   Data: df0
Control: control

REML criterion at convergence: 24351.7

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1105 -0.5636 -0.0071  0.5640  7.3716 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Name     (Intercept) 0.006226 0.07891      
          Time        0.004116 0.06415  0.18
 Residual             0.066142 0.25718      
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                   Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)        0.505972      0.004329   4392.427008 116.881 < 0.0000000000000002 ***
Time              -0.148762      0.016674 151068.334259  -8.922 < 0.0000000000000002 ***
PeriodPB           0.022615      0.003709 169246.914124   6.098        0.00000000108 ***
PeriodPA           0.022052      0.005557 169331.953373   3.968        0.00007248943 ***
Time:PeriodPB      0.193346      0.016754 169082.738995  11.540 < 0.0000000000000002 ***
Time:PeriodPA      0.083065      0.017673 168894.566947   4.700        0.00000260350 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Time   PerdPB PerdPA Tm:PPB
Time        -0.635                            
PeriodPB    -0.699  0.759                     
PeriodPA    -0.461  0.512  0.536              
Time:PerdPB  0.647 -0.976 -0.692 -0.509       
Time:PerdPA  0.616 -0.923 -0.718 -0.759  0.918

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

# R2 for Mixed Models

  Conditional R2: 0.105
     Marginal R2: 0.004

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.101
  Unadjusted ICC: 0.101

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

Model contains random slopes. Cannot compute accurate ICCs by group factors.

# ICC by Group

Group |   ICC
-------------
Name  | 0.085

Model Summary

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

NOTE: Results may be misleading due to involvement in interactions

tbl0

Characteristic	Beta	95% CI1	p-value
Time	-0.15	-0.18, -0.12	<0.001
Period
PB - PC	0.01	0.00, 0.02	0.10
PA - PC	0.02	0.00, 0.03	0.032
PA - PB	0.01	0.00, 0.02	0.4
Time * Period
Time * PB	0.19	0.16, 0.23	<0.001
Time * PA	0.08	0.05, 0.12	<0.001
Name.sd__(Intercept)	0.08
Name.cor__(Intercept).Time	0.18
Name.sd__Time	0.06
Residual.sd__Observation	0.26
1 CI = Confidence Interval

Model Parameters

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

parameters::model_parameters(get(model))

# Fixed Effects

Parameter          | Coefficient |       SE |         95% CI | t(169987) |      p
---------------------------------------------------------------------------------
(Intercept)        |        0.51 | 4.33e-03 | [ 0.50,  0.51] |    116.88 | < .001
Time               |       -0.15 |     0.02 | [-0.18, -0.12] |     -8.92 | < .001
Period [PB]        |        0.02 | 3.71e-03 | [ 0.02,  0.03] |      6.10 | < .001
Period [PA]        |        0.02 | 5.56e-03 | [ 0.01,  0.03] |      3.97 | < .001
Time × Period [PB] |        0.19 |     0.02 | [ 0.16,  0.23] |     11.54 | < .001
Time × Period [PA] |        0.08 |     0.02 | [ 0.05,  0.12] |      4.70 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.08
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |        0.18
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

Plot Model: Estimates

pp(model)

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

extra <- "0000"
terms <- "NONE"

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(get(model), type="est")
gg88 <- gg88 + ggplot2::ylim(-0.3,0.3) + line2 + stuff0

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-estimates-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Period

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

extra <- "0002"
terms=c("Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Period

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

extra <- "0003"
terms=c("Time", "Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time

NULL: Time

pp(model)

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.53, 0.55
-0.50 |      0.53 | 0.52, 0.54
 0.00 |      0.53 | 0.52, 0.53
 0.50 |      0.52 | 0.52, 0.53
 1.00 |      0.52 | 0.51, 0.53

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Slope     |       95% CI |     p
--------------------------------
-9.10e-03 | -0.02,  0.00 | 0.025

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,) 
gg88 <- gg88 + line0 + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Period

NULL: Period

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

extra <- "1002"
terms <- c("Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Period | Predicted |     95% CI
-------------------------------
PC     |      0.52 | 0.51, 0.53
PB     |      0.53 | 0.52, 0.53
PA     |      0.54 | 0.52, 0.55

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Period | Predicted |     95% CI |      p
----------------------------------------
PB     |      0.53 | 0.52, 0.53 | < .001
PC     |      0.52 | 0.51, 0.53 | < .001
PA     |      0.54 | 0.52, 0.55 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Period

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
## gg88 <- gg88 + line0 + line1 + time0 + stuff0 

gg88 <- gg88 + stuff0 

gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Period

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Period |  Contrast |       95% CI |     p
-----------------------------------------
PB-PC  |  9.46e-03 |  0.00,  0.02 | 0.058
PB-PA  | -6.94e-03 | -0.02,  0.00 | 0.173
PC-PA  |     -0.02 | -0.03,  0.00 | 0.036

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Period

NULL: Time x Period

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

extra <- "1003"
terms <- c("Time", "Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")


cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.65 | 0.61, 0.69
-0.50 |      0.58 | 0.56, 0.60
 0.00 |      0.51 | 0.50, 0.52
 0.50 |      0.44 | 0.42, 0.45
 1.00 |      0.36 | 0.34, 0.39

Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.47, 0.49
-0.50 |      0.51 | 0.50, 0.51
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.56 | 0.55, 0.56
 1.00 |      0.58 | 0.57, 0.59

Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.59 | 0.57, 0.61
-0.50 |      0.56 | 0.55, 0.58
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.50 | 0.49, 0.51
 1.00 |      0.47 | 0.46, 0.48

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Period | Slope |       95% CI |      p
--------------------------------------
PC     | -0.15 | -0.18, -0.12 | < .001
PB     |  0.04 |  0.04,  0.05 | < .001
PA     | -0.07 | -0.08, -0.05 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 
gg88 <- gg88 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Period

pp(model);

fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Period | Contrast |       95% CI |      p
-----------------------------------------
PC-PB  |    -0.19 | -0.23, -0.16 | < .001
PC-PA  |    -0.08 | -0.12, -0.05 | < .001
PB-PA  |     0.11 |  0.10,  0.12 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

Model 04: Time x Period x Outcome

Fit

model <- "fit04aPd"
suppressWarnings(rm(list = model))
assign(
  model,
  lmerTest::lmer(
    formula = Agency ~ (Time | Name) + Time * Period * Outcome,
    data = df0,
    REML = REML,
    control = control))

fbase <- file.path(dir8, model); write_model_info(fbase)

readr::write_rds(get(model), file = paste0(fbase, ".mdl.rds"))
pp(model); summary(get(model)) %>% print(correlation = TRUE)

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Agency ~ (Time | Name) + Time * Period * Outcome
   Data: df0
Control: control

REML criterion at convergence: 23854.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-8.1837 -0.5638 -0.0081  0.5631  7.3819 

Random effects:
 Groups   Name        Variance Std.Dev. Corr 
 Name     (Intercept) 0.004807 0.06934       
          Time        0.003341 0.05780  -0.07
 Residual             0.066045 0.25699       
Number of obs: 169997, groups:  Name, 870

Fixed effects:
                                 Estimate    Std. Error            df t value             Pr(>|t|)    
(Intercept)                      0.456679      0.007401  10611.863226  61.707 < 0.0000000000000002 ***
Time                            -0.270178      0.032339 160626.153738  -8.355 < 0.0000000000000002 ***
PeriodPB                         0.058436      0.006870 169298.359477   8.506 < 0.0000000000000002 ***
PeriodPA                        -0.046169      0.011403 168642.624739  -4.049    0.000051518894422 ***
Outcomewinner                    0.077669      0.008954   7522.582122   8.674 < 0.0000000000000002 ***
Time:PeriodPB                    0.327201      0.032499 169002.252327  10.068 < 0.0000000000000002 ***
Time:PeriodPA                    0.250601      0.034945 169220.044400   7.171    0.000000000000746 ***
Time:Outcomewinner               0.157067      0.037712 158776.211635   4.165    0.000031164035167 ***
PeriodPB:Outcomewinner          -0.053725      0.008169 169358.119965  -6.577    0.000000000048304 ***
PeriodPA:Outcomewinner           0.086802      0.013062 169180.936257   6.645    0.000000000030345 ***
Time:PeriodPB:Outcomewinner     -0.177453      0.037926 169245.190328  -4.679    0.000002885388234 ***
Time:PeriodPA:Outcomewinner     -0.213600      0.040515 169254.780591  -5.272    0.000000135027103 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) Time   PerdPB PerdPA Otcmwn Tm:PPB Tm:PPA Tm:Otc PrPB:O PrPA:O T:PPB:
Time        -0.715                                                                      
PeriodPB    -0.808  0.769                                                               
PeriodPA    -0.478  0.471  0.514                                                        
Outcomewnnr -0.827  0.591  0.668  0.395                                                 
Time:PerdPB  0.709 -0.985 -0.722 -0.469 -0.586                                          
Time:PerdPA  0.662 -0.912 -0.713 -0.750 -0.547  0.907                                   
Tm:Otcmwnnr  0.613 -0.858 -0.659 -0.404 -0.701  0.845  0.782                            
PrdPB:Otcmw  0.679 -0.646 -0.841 -0.432 -0.770  0.607  0.600  0.765                     
PrdPA:Otcmw  0.417 -0.411 -0.448 -0.873 -0.475  0.410  0.655  0.485  0.519              
Tm:PrdPB:Ot -0.608  0.844  0.619  0.402  0.693 -0.857 -0.778 -0.982 -0.708 -0.482       
Tm:PrdPA:Ot -0.571  0.786  0.615  0.647  0.651 -0.783 -0.863 -0.916 -0.713 -0.753  0.911

cat0(sep2); pp(model); performance::r2(get(model))

--------------------------------------------------------------------- 

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

# R2 for Mixed Models

  Conditional R2: 0.100
     Marginal R2: 0.019

cat0(sep2); pp(model); performance::icc(get(model))

--------------------------------------------------------------------- 

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.083
  Unadjusted ICC: 0.081

cat0(sep2); pp(model); performance::icc(get(model), by_group = TRUE)

--------------------------------------------------------------------- 

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

Model contains random slopes. Cannot compute accurate ICCs by group factors.

# ICC by Group

Group |   ICC
-------------
Name  | 0.067

Model Summary

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

tbl0 <- gtsummary::tbl_regression(
  get(model),
  add_pairwise_contrasts = TRUE,
  exponentiate = FALSE,
  tidy_fun = broom.mixed::tidy,
)

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

NOTE: Results may be misleading due to involvement in interactions

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'pbkrtest.limit = 169997' (or larger)
[or, globally, 'set emm_options(pbkrtest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

Note: D.f. calculations have been disabled because the number of observations exceeds 3000.
To enable adjustments, add the argument 'lmerTest.limit = 169997' (or larger)
[or, globally, 'set emm_options(lmerTest.limit = 169997)' or larger];
but be warned that this may result in large computation time and memory use.

NOTE: Results may be misleading due to involvement in interactions

tbl0

Characteristic	Beta	95% CI1	p-value
Time	-0.27	-0.33, -0.21	<0.001
Period
PB - PC	0.02	0.00, 0.03	0.007
PA - PC	-0.01	-0.03, 0.01	0.2
PA - PB	-0.03	-0.04, -0.01	<0.001
Outcome
winner - loser	0.09	0.07, 0.10	<0.001
Time * Period
Time * PB	0.33	0.26, 0.39	<0.001
Time * PA	0.25	0.18, 0.32	<0.001
Time * Outcome
Time * winner	0.16	0.08, 0.23	<0.001
Period * Outcome
PB * winner	-0.05	-0.07, -0.04	<0.001
PA * winner	0.09	0.06, 0.11	<0.001
Time * Period * Outcome
Time * PB * winner	-0.18	-0.25, -0.10	<0.001
Time * PA * winner	-0.21	-0.29, -0.13	<0.001
Name.sd__(Intercept)	0.07
Name.cor__(Intercept).Time	-0.07
Name.sd__Time	0.06
Residual.sd__Observation	0.26
1 CI = Confidence Interval

Model Parameters

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

parameters::model_parameters(get(model))

# Fixed Effects

Parameter                               | Coefficient |       SE |         95% CI | t(169981) |      p
------------------------------------------------------------------------------------------------------
(Intercept)                             |        0.46 | 7.40e-03 | [ 0.44,  0.47] |     61.71 | < .001
Time                                    |       -0.27 |     0.03 | [-0.33, -0.21] |     -8.35 | < .001
Period [PB]                             |        0.06 | 6.87e-03 | [ 0.04,  0.07] |      8.51 | < .001
Period [PA]                             |       -0.05 |     0.01 | [-0.07, -0.02] |     -4.05 | < .001
Outcome [winner]                        |        0.08 | 8.95e-03 | [ 0.06,  0.10] |      8.67 | < .001
Time × Period [PB]                      |        0.33 |     0.03 | [ 0.26,  0.39] |     10.07 | < .001
Time × Period [PA]                      |        0.25 |     0.03 | [ 0.18,  0.32] |      7.17 | < .001
Time × Outcome [winner]                 |        0.16 |     0.04 | [ 0.08,  0.23] |      4.16 | < .001
Period [PB] × Outcome [winner]          |       -0.05 | 8.17e-03 | [-0.07, -0.04] |     -6.58 | < .001
Period [PA] × Outcome [winner]          |        0.09 |     0.01 | [ 0.06,  0.11] |      6.65 | < .001
(Time × Period [PB]) × Outcome [winner] |       -0.18 |     0.04 | [-0.25, -0.10] |     -4.68 | < .001
(Time × Period [PA]) × Outcome [winner] |       -0.21 |     0.04 | [-0.29, -0.13] |     -5.27 | < .001

# Random Effects

Parameter                  | Coefficient
----------------------------------------
SD (Intercept: Name)       |        0.07
SD (Time: Name)            |        0.06
Cor (Intercept~Time: Name) |       -0.07
SD (Residual)              |        0.26

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed using a Wald t-distribution approximation.

FE: Time

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "0001"
terms=c("Time")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + time0 + stuff0 

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Period

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "0002"
terms=c("Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Period

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "0003"
terms=c("Time", "Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "0004"
terms=c("Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "0005"
terms=c("Time", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + line0 + line1 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Period x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "0006"
terms=c("Period", "Outcome")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88 <- gg88 + stuff0

gg88$data$x <- c(2,1,3) ## CAUTION: THIS WORKS WITH: sjPlot (OR SIMPLER MODEL)
## gg88$data$x <- c(2,2,1,1,3,3) ## CAUTION: THIS WORKS WITH: ggeffects
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

FE: Time x Outcome x Period

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "0007"
terms=c("Time", "Outcome", "Period")

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
  get(model), terms=terms, type="eff", show.data=FALSE) 

gg88$data$facet <- factor(
  case_when(
    gg88$data$facet == "PB" ~ "PB",
    gg88$data$facet == "PC" ~ "PC",
    gg88$data$facet == "PA" ~ "PA",
  ),
  levels=c("PB", "PC", "PA"))

gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect0

## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)
## gg88 <- gg88 + ggplot2::coord_cartesian(ylim = c(0.25, 0.75))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".png"))
ggsave(file=file, plot=gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-fixed-effects-", paste(terms, collapse="-x-"),".svg"))
ggsave(file=file, plot=gg88, width=8, height=6)
gg88

EF: Time:

NULL: Time

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1001"
terms <- c("Time")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.54 | 0.53, 0.55
-0.50 |      0.53 | 0.53, 0.54
 0.00 |      0.53 | 0.52, 0.53
 0.50 |      0.52 | 0.52, 0.53
 1.00 |      0.52 | 0.51, 0.53

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Slope |       95% CI |     p
----------------------------
-0.01 | -0.02,  0.00 | 0.005

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05,)

gg88 <- gg88 + line0 + line1 + time0 + stuff0 
## gg88 <- gg88 + coord_cartesian(ylim = c(-0.25, 1.25))

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

EF: Period

NULL: Period

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1002"
terms <- c("Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Period | Predicted |     95% CI
-------------------------------
PC     |      0.52 | 0.51, 0.54
PB     |      0.52 | 0.52, 0.53
PA     |      0.52 | 0.51, 0.53

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Period | Predicted |     95% CI |      p
----------------------------------------
PB     |      0.52 | 0.52, 0.53 | < .001
PC     |      0.52 | 0.51, 0.54 | < .001
PA     |      0.52 | 0.51, 0.53 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Period

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(
  limit_range=FALSE,  
  show_data = FALSE, 
  dot_alpha = 0.05)

gg88 <- gg88 + stuff0
gg88$data$x <- c(2,1,3)
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Period

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Period | Contrast |      95% CI |     p
---------------------------------------
PB-PC  | 3.02e-04 | -0.01, 0.01 | 0.952
PB-PA  | 4.20e-03 | -0.01, 0.01 | 0.874
PC-PA  | 3.90e-03 | -0.01, 0.02 | 0.874

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Period

NULL: Time x Period

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1003"
terms <- c("Period", "Time")
terms <- c("Time", "Period")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.67 | 0.63, 0.71
-0.50 |      0.59 | 0.56, 0.61
 0.00 |      0.51 | 0.50, 0.51
 0.50 |      0.43 | 0.41, 0.44
 1.00 |      0.34 | 0.32, 0.37

Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.48, 0.49
-0.50 |      0.51 | 0.50, 0.51
 0.00 |      0.53 | 0.52, 0.53
 0.50 |      0.55 | 0.54, 0.56
 1.00 |      0.57 | 0.56, 0.58

Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.57 | 0.55, 0.60
-0.50 |      0.55 | 0.53, 0.56
 0.00 |      0.52 | 0.51, 0.53
 0.50 |      0.49 | 0.49, 0.50
 1.00 |      0.46 | 0.46, 0.47

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Period | Slope |       95% CI |      p
--------------------------------------
PC     | -0.27 | -0.33, -0.21 | < .001
PB     |  0.06 |  0.05,  0.07 | < .001
PA     | -0.02 | -0.05,  0.01 | 0.173 

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect2

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Time x Period

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Period | Contrast |       95% CI |      p
-----------------------------------------
PC-PB  |    -0.33 | -0.39, -0.26 | < .001
PC-PA  |    -0.25 | -0.32, -0.18 | < .001
PB-PA  |     0.08 |  0.05,  0.11 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Outcome

NULL: Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1004"
terms <- c("Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome | Predicted |     95% CI
--------------------------------
loser   |      0.45 | 0.44, 0.46
winner  |      0.52 | 0.51, 0.52

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Outcome | Predicted |     95% CI |      p
-----------------------------------------
loser   |      0.45 | 0.44, 0.46 | < .001
winner  |      0.52 | 0.51, 0.52 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Outcome

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
ggsave(file = "file.svg", plot = gg88, width=8, height=6)
gg88

DIFF: Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Outcome      | Contrast |       95% CI |      p
-----------------------------------------------
loser-winner |    -0.07 | -0.08, -0.06 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Outcome:

NULL: Time x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1005"
terms <- c("Time", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome: loser

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.48 | 0.46, 0.50
-0.50 |      0.48 | 0.47, 0.49
 0.00 |      0.47 | 0.47, 0.48
 0.50 |      0.47 | 0.46, 0.48
 1.00 |      0.47 | 0.45, 0.48

Outcome: winner

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.56 | 0.55, 0.57
-0.50 |      0.55 | 0.54, 0.56
 0.00 |      0.55 | 0.54, 0.55
 0.50 |      0.54 | 0.53, 0.55
 1.00 |      0.53 | 0.52, 0.54

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Outcome | Slope |       95% CI |      p
---------------------------------------
loser   | -0.27 | -0.33, -0.21 | < .001
winner  | -0.11 | -0.15, -0.08 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT Time x Outcome

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)
gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect3

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

DIFF: Time x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Outcome      | Contrast |       95% CI |      p
-----------------------------------------------
loser-winner |    -0.16 | -0.23, -0.08 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Period x Outcome

NULL: Period x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1006"
terms <- c("Period", "Outcome")
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome: loser

Period | Predicted |     95% CI
-------------------------------
PC     |      0.47 | 0.45, 0.49
PB     |      0.51 | 0.50, 0.52
PA     |      0.41 | 0.39, 0.43

Outcome: winner

Period | Predicted |     95% CI
-------------------------------
PC     |      0.54 | 0.53, 0.55
PB     |      0.53 | 0.53, 0.54
PA     |      0.58 | 0.56, 0.59

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

Period | Outcome | Predicted |     95% CI |      p
--------------------------------------------------
PC     |   loser |      0.47 | 0.45, 0.49 | < .001
PB     |   loser |      0.51 | 0.50, 0.52 | < .001
PA     |   loser |      0.41 | 0.39, 0.43 | < .001
PC     |  winner |      0.54 | 0.53, 0.55 | < .001
PB     |  winner |      0.53 | 0.53, 0.54 | < .001
PA     |  winner |      0.58 | 0.56, 0.59 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Period x Outcome

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

## gg88$data$x <- c(2,1,3) ## CAUTION: THIS WORKS WITH sjPlot (OR SIMPLER MODEL)
gg88$data$x <- c(2,2,1,1,3,3) ## CAUTION: THIS WORKS WITH: ggeffects
gg88 <- gg88 + scale_x_continuous(breaks=c(1,2,3), labels=c("PB", "PC", "PA"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

gg88 <- gg88 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=8, height=6)
gg88

PLOT: Period x Outcome [low-level]

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1006"
terms <- c("Period", "Outcome")

type <- "eff"
type <- "pred"

suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- sjPlot::plot_model(
  get(model), 
  type = type, 
  terms = terms, 
  )

## str(gg44$data)
gg44$data$x <- c(2,2,1,1,3,3)
## gg44 <- gg44 + scale_x_discrete(labels = c("BE", "MI", "ARRR"))
gg44 <- gg44 + scale_x_continuous(
  breaks=c(1,2,3),
  labels=c("PB", "PC", "PA"))

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

gg44 <- gg44 + stuff0

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg44, width=8, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjplot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg44, width=8, height=6)
gg44

DIFF: Period x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# Pairwise comparisons

Period |       Outcome | Contrast |       95% CI |      p
---------------------------------------------------------
PC-PB  |   loser-loser |    -0.04 | -0.05, -0.02 | < .001
PC-PA  |   loser-loser |     0.06 |  0.04,  0.09 | < .001
PC-PC  |  loser-winner |    -0.07 | -0.09, -0.05 | < .001
PC-PB  |  loser-winner |    -0.06 | -0.08, -0.04 | < .001
PC-PA  |  loser-winner |    -0.11 | -0.13, -0.08 | < .001
PB-PA  |   loser-loser |     0.10 |  0.08,  0.12 | < .001
PB-PC  |  loser-winner |    -0.03 | -0.05, -0.02 | < .001
PB-PB  |  loser-winner |    -0.03 | -0.04, -0.01 | < .001
PB-PA  |  loser-winner |    -0.07 | -0.08, -0.05 | < .001
PA-PC  |  loser-winner |    -0.13 | -0.15, -0.11 | < .001
PA-PB  |  loser-winner |    -0.12 | -0.15, -0.10 | < .001
PA-PA  |  loser-winner |    -0.17 | -0.19, -0.14 | < .001
PC-PB  | winner-winner | 5.48e-03 | -0.01,  0.02 | 0.311 
PC-PA  | winner-winner |    -0.04 | -0.05, -0.02 | < .001
PB-PA  | winner-winner |    -0.04 | -0.06, -0.03 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

EF: Time x Period x Outcome

NULL: Time x Period x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

extra <- "1007"
terms <- c("Time", "Period", "Outcome")
terms <- c("Time", "Outcome", "Period") ## CAUTION
pred0 <- ggeffects::predict_response(get(model), terms=terms)
test0 <- ggeffects::test_predictions(get(model), terms=terms, test=NULL, p_adjust="fdr")

cat0(sep0)

===================================================================== 

pred0 %>% print(n=Inf)

# Average predicted values of Agency

Outcome: loser
Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.72 | 0.65, 0.80
-0.50 |      0.59 | 0.54, 0.63
 0.00 |      0.45 | 0.44, 0.47
 0.50 |      0.32 | 0.30, 0.34
 1.00 |      0.19 | 0.13, 0.24

Outcome: loser
Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.45 | 0.44, 0.46
-0.50 |      0.48 | 0.47, 0.49
 0.00 |      0.51 | 0.50, 0.52
 0.50 |      0.54 | 0.53, 0.55
 1.00 |      0.57 | 0.56, 0.59

Outcome: loser
Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.43 | 0.38, 0.47
-0.50 |      0.42 | 0.38, 0.45
 0.00 |      0.41 | 0.39, 0.43
 0.50 |      0.40 | 0.39, 0.41
 1.00 |      0.39 | 0.38, 0.41

Outcome: winner
Period: PC

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.64 | 0.60, 0.69
-0.50 |      0.59 | 0.56, 0.61
 0.00 |      0.53 | 0.52, 0.54
 0.50 |      0.48 | 0.46, 0.49
 1.00 |      0.42 | 0.39, 0.45

Outcome: winner
Period: PB

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.50 | 0.49, 0.51
-0.50 |      0.52 | 0.51, 0.52
 0.00 |      0.54 | 0.53, 0.54
 0.50 |      0.56 | 0.55, 0.57
 1.00 |      0.58 | 0.56, 0.59

Outcome: winner
Period: PA

 Time | Predicted |     95% CI
------------------------------
-1.00 |      0.65 | 0.62, 0.67
-0.50 |      0.61 | 0.59, 0.63
 0.00 |      0.57 | 0.56, 0.58
 0.50 |      0.54 | 0.53, 0.54
 1.00 |      0.50 | 0.49, 0.51

cat0(sep0)

===================================================================== 

test0 %>% print(n=Inf)

# (Average) Linear trend for Time

Outcome | Period | Slope |       95% CI |      p
------------------------------------------------
loser   |     PC | -0.27 | -0.33, -0.21 | < .001
loser   |     PB |  0.06 |  0.05,  0.07 | < .001
loser   |     PA | -0.02 | -0.05,  0.01 | 0.173 
winner  |     PC | -0.11 | -0.15, -0.08 | < .001
winner  |     PB |  0.04 |  0.03,  0.05 | < .001
winner  |     PA | -0.08 | -0.09, -0.06 | < .001

pred0tab <- ggeffects::print_html(pred0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.html"))
pred0tab %>% tinytable::save_tt(file, overwrite=TRUE)

test0tab <- ggeffects::print_html(test0)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test0.html"))
test0tab %>% tinytable::save_tt(file, overwrite=TRUE)

PLOT: Time x Period x Outcome

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- pred0 %>% plot(limit_range=FALSE,  show_data = FALSE, dot_alpha = 0.05)

gg88$data$facet <- factor(
  case_when(
    gg88$data$facet == "PB" ~ "PB",
    gg88$data$facet == "PC" ~ "PC",
    gg88$data$facet == "PA" ~ "PA",
  ),
  levels=c("PB", "PC", "PA"))

gg88 <- gg88 + line0 + line1 + time0 + stuff0 

gg88 <- gg88 + rect0

## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)

file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=12, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=12, height=6)
gg88

PLOT: Time x Period x Outcome (low-level)

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_model(
          get(model),
          type="emm",
          terms=terms) 

gg88$data$facet <- factor(
  case_when(
    gg88$data$facet == "PB" ~ "PB",
    gg88$data$facet == "PC" ~ "PC",
    gg88$data$facet == "PA" ~ "PA",
  ),
  levels=c("PB", "PC", "PA"))

gg88 <- gg88 + line0 + line1 + time0 + stuff0 + rect0
## gg88 <- gg88 + ggplot2::ylim(0.15,0.85)

file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.png"))
ggsave(file = file, plot = gg88, width=12, height=6)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-sjPlot-", paste(terms, collapse="-x-"),"-pred0.svg"))
ggsave(file = file, plot = gg88, width=12, height=6)
gg88

DIFF: Time x Period x Outcome

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

test2 <- ggeffects::test_predictions(get(model), terms=terms, test="pairwise", p_adjust="fdr")
test2 %>% print()

# (Average) Linear trend for Time

Outcome       | Period | Contrast |       95% CI |      p
---------------------------------------------------------
loser-loser   |  PC-PB |    -0.33 | -0.39, -0.26 | < .001
loser-loser   |  PC-PA |    -0.25 | -0.32, -0.18 | < .001
loser-winner  |  PC-PC |    -0.16 | -0.23, -0.08 | < .001
loser-winner  |  PC-PB |    -0.31 | -0.37, -0.24 | < .001
loser-winner  |  PC-PA |    -0.19 | -0.26, -0.13 | < .001
loser-loser   |  PB-PA |     0.08 |  0.05,  0.11 | < .001
loser-winner  |  PB-PC |     0.17 |  0.13,  0.21 | < .001
loser-winner  |  PB-PB |     0.02 |  0.01,  0.03 | 0.005 
loser-winner  |  PB-PA |     0.13 |  0.11,  0.15 | < .001
loser-winner  |  PA-PC |     0.09 |  0.05,  0.14 | < .001
loser-winner  |  PA-PB |    -0.06 | -0.09, -0.03 | < .001
loser-winner  |  PA-PA |     0.06 |  0.02,  0.09 | < .001
winner-winner |  PC-PB |    -0.15 | -0.19, -0.11 | < .001
winner-winner |  PC-PA |    -0.04 | -0.08,  0.00 | 0.071 
winner-winner |  PB-PA |     0.11 |  0.10,  0.13 | < .001

test2tab <- ggeffects::print_html(test2)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.html"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)
file = file.path(dir8, paste0(model,"-xtr-",extra,"-ggeffects-", paste(terms, collapse="-x-"),"-test2.docx"))
test2tab %>% tinytable::save_tt(file, overwrite=TRUE)

CHECK: Model Assumptions

pp(model);

fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

performance::check_model(get(model))

ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)

Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.

ggsave(
  file = "performance-test.png", 
  ## file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  ## plot = gg88,
  ## plot = last_plot(),
  width=8,
  height=12)

Warning: Failed to fit group -1.
Failed to fit group -1.
Caused by error in `predLoess()`:
! workspace required (43350042532) is too large probably because of setting 'se = TRUE'.

ls()

 [1] "cat0"                     "check_pkgs"               "control"                  "create_factor_mappings"   "df0"                      "df2"                      "df3"                      "dir8"                    
 [9] "extra"                    "fbase"                    "file"                     "fit01aPd"                 "fit02aPd"                 "fit03aPd"                 "fit04aPd"                 "fpath"                   
[17] "generate_permutations"    "get_stars"                "gg44"                     "gg88"                     "ggpairs_lower_fun"        "hasConverged"             "ifd0"                     "ifn0"                    
[25] "jiko_count_na"            "line0"                    "line1"                    "line2"                    "model"                    "name"                     "obj_has_name"             "ofd0"                    
[33] "ofd2"                     "pp"                       "pred0"                    "pred0tab"                 "rect0"                    "rect2"                    "rect3"                    "REML"                    
[41] "replace_factor_levels"    "replace_outliers_with_NA" "rescale01"                "rxy0"                     "rxy2"                     "rxy3"                     "sep0"                     "sep1"                    
[49] "sep2"                     "set.w"                    "stuff0"                   "tbl0"                     "terms"                    "test0"                    "test0tab"                 "test2"                   
[57] "test2tab"                 "time0"                    "type"                     "var0"                     "write_model_info"        

Auxiliary Checkups

Variables influence / importance in fit04aPd

## TODO CHANGE TITLE
formula(fit04aPd)

Agency ~ (Time | Name) + Time * Period * Outcome

eff_size0 <- effectsize::eta_squared(fit04aPd)
eff_size2 <- eff_size0 %>%
  dplyr::mutate(Interpret = effectsize::interpret_eta_squared(Eta2_partial))

performance::print_html(eff_size2)

Effect Size for ANOVA
Parameter	Eta2 (partial)	95% CI	Interpret
Time	2.12e-03	[0.00, 1.00]	very small
Period	5.13e-04	[0.00, 1.00]	very small
Outcome	0.07	[0.05, 1.00]	medium
Time:Period	1.65e-03	[0.00, 1.00]	very small
Time:Outcome	9.19e-05	[0.00, 1.00]	very small
Period:Outcome	1.08e-03	[0.00, 1.00]	very small
Time:Period:Outcome	1.71e-04	[0.00, 1.00]	very small

List models (for Model Hierarchy / Hierarchies)

temp <- lapply(ls(pattern="^fit\\d+"), pp)

fit01aPd: [df0] Agency ~ (1 | Name) + 1
fit02aPd: [df0] Agency ~ (Time | Name) + Time
fit03aPd: [df0] Agency ~ (Time | Name) + Time * Period
fit04aPd: [df0] Agency ~ (Time | Name) + Time * Period * Outcome

Compare models using performance package with score

perf0 <- performance::compare_performance(
  fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf0 %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit04aPd	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	84.36%
fit03aPd	lmerModLmerTest	0.11	4.32e-03	0.10	0.26	0.26	4.74e-116	4.74e-116	5.77e-103	51.23%
fit02aPd	lmerModLmerTest	0.10	4.62e-04	0.10	0.26	0.26	5.28e-250	5.29e-250	3.40e-228	43.11%
fit01aPd	lmerModLmerTest	0.08	0.00	0.08	0.26	0.26	0.00e+00	0.00e+00	0.00e+00	1.03%
NA

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- perf0 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf0.compare-02-score.png"),
  plot = gg88,
  width=8,
  height=6)

gg88

Compare less models

perf2 <- performance::compare_performance(
  ## fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  ## fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
  ## CAUTION: COMMA
  rank = TRUE, verbose = FALSE)
perf2 %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit04aPd	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	75.00%
fit03aPd	lmerModLmerTest	0.11	4.32e-03	0.10	0.26	0.26	4.74e-116	4.74e-116	5.77e-103	25.00%
NA

suppressWarnings(rm(list = ls(pattern = "^gg44")))
gg44 <- perf2 %>% plot()
ggsave(
  file = file.path(dir8, "summary-combined-models.perf2.compare-03-score.png"),
  plot = gg44,
  width=8,
  height=6)

gg44  

Comment

Add interpretation of the difference between best fit to data and best predictive power and overfitting (aspecially in the context of the last model).

gg44

Performance Table Sorted by R2_conditional

perf0 %>% dplyr::arrange(desc(R2_conditional)) %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit03aPd	lmerModLmerTest	0.11	4.32e-03	0.10	0.26	0.26	4.74e-116	4.74e-116	5.77e-103	51.23%
fit02aPd	lmerModLmerTest	0.10	4.62e-04	0.10	0.26	0.26	5.28e-250	5.29e-250	3.40e-228	43.11%
fit04aPd	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	84.36%
fit01aPd	lmerModLmerTest	0.08	0.00	0.08	0.26	0.26	0.00e+00	0.00e+00	0.00e+00	1.03%
NA

Performance Table Sorted by R2_marginal

perf0 %>% dplyr::arrange(desc(R2_marginal)) %>% performance::print_html()

Comparison of Model Performance Indices
Name	Model	R2 (cond.)	R2 (marg.)	ICC	RMSE	Sigma	AIC weights	AICc weights	BIC weights	Performance-Score
fit04aPd	lmerModLmerTest	0.10	0.02	0.08	0.26	0.26	1.00	1.00	1.00	84.36%
fit03aPd	lmerModLmerTest	0.11	4.32e-03	0.10	0.26	0.26	4.74e-116	4.74e-116	5.77e-103	51.23%
fit02aPd	lmerModLmerTest	0.10	4.62e-04	0.10	0.26	0.26	5.28e-250	5.29e-250	3.40e-228	43.11%
fit01aPd	lmerModLmerTest	0.08	0.00	0.08	0.26	0.26	0.00e+00	0.00e+00	0.00e+00	1.03%
NA

Interpret R2

model <- "fit01aPd" # [df0] Agency ~ (1 | Name) + 1
model <- "fit02aPd" # [df0] Agency ~ (Time | Name) + Time
model <- "fit03aPd" # [df0] Agency ~ (Time | Name) + Time * Period
model <- "fit04aPd" # [df0] Agency ~ (Time | Name) + Time * Period * Outcome

cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_conditional, rules="cohen1988"))

weak 

cat0(effectsize::interpret_r2(performance::r2(get(model))$R2_marginal, rules="cohen1988"))

very weak 

Plot models: Part 1: Basic Models

suppressWarnings(rm(list = ls(pattern = "^gg88")))
gg88 <- sjPlot::plot_models(
  ## CAUTION the null model can not be used here 
  ## Thus to keep the numbers consistent I have 
  ## used model 02 as an input twice
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time 
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
  spacing=1, dot.size=1
)

ggsave(
  file = file.path(dir8, "summary-all-models-part-01.png"),
  plot = gg88,
  width=5,
  height=4)

gg88

Tab Models

The above `plot_models` seems to be a lot more illustrative

# , eval = TRUE, results='hide'
library(sjPlot)
library(sjmisc)

Attaching package: 'sjmisc'

The following object is masked from 'package:skimr':

    to_long

The following object is masked from 'package:purrr':

    is_empty

The following object is masked from 'package:tidyr':

    replace_na

The following object is masked from 'package:tibble':

    add_case

library(sjlabelled)

Attaching package: 'sjlabelled'

The following object is masked from 'package:forcats':

    as_factor

The following object is masked from 'package:dplyr':

    as_label

The following object is masked from 'package:ggplot2':

    as_label

tab0 <- sjPlot::tab_model(
  fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
    show.reflvl = TRUE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-01-SEE-PNG-images-for-a-better-overview.html"))
tab0

	Agency			Agency			Agency			Agency
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	0.50 ***	0.49 – 0.50	<0.001	0.50 ***	0.49 – 0.50	<0.001	0.51 ***	0.50 – 0.51	<0.001	0.46 ***	0.44 – 0.47	<0.001
PC	Reference			Reference			Reference			Reference
PB							0.02 ***	0.02 – 0.03	<0.001	0.06 ***	0.04 – 0.07	<0.001
PA							0.02 ***	0.01 – 0.03	<0.001	-0.05 ***	-0.07 – -0.02	<0.001
PeriodPA:Outcomewinner										0.09 ***	0.06 – 0.11	<0.001
loser	Reference			Reference			Reference			Reference
PeriodPB:Outcomewinner										-0.05 ***	-0.07 – -0.04	<0.001
winner										0.08 ***	0.06 – 0.10	<0.001
Time				-0.01 ***	-0.02 – -0.01	<0.001	-0.15 ***	-0.18 – -0.12	<0.001	-0.27 ***	-0.33 – -0.21	<0.001
Time:Outcomewinner										0.16 ***	0.08 – 0.23	<0.001
Time:PeriodPA							0.08 ***	0.05 – 0.12	<0.001	0.25 ***	0.18 – 0.32	<0.001
Time:PeriodPA:Outcomewinner										-0.21 ***	-0.29 – -0.13	<0.001
Time:PeriodPB							0.19 ***	0.16 – 0.23	<0.001	0.33 ***	0.26 – 0.39	<0.001
Time:PeriodPB:Outcomewinner										-0.18 ***	-0.25 – -0.10	<0.001
Random Effects
σ2	0.07			0.07			0.07			0.07
τ00	0.01 Name			0.01 Name			0.01 Name			0.00 Name
τ11				0.00 Name.Time			0.00 Name.Time			0.00 Name.Time
ρ01				0.18 Name			0.18 Name			-0.07 Name
ICC	0.08			0.10			0.10			0.08
N	870 Name			870 Name			870 Name			870 Name
Observations	169997			169997			169997			169997
Marginal R2 / Conditional R2	0.000 / 0.084			0.000 / 0.102			0.004 / 0.105			0.019 / 0.100
* p<0.05   ** p<0.01   *** p<0.001

## knitr::html
## papaja::
## rmarkdown::render
## rmarkdown::render(tab0$knitr)
## tab0$knitr

Tabulate models (publication ready)

tab0 <- sjPlot::tab_model(
  fit01aPd, # [df0] Agency ~ (1 | Name) + 1
  fit02aPd, # [df0] Agency ~ (Time | Name) + Time
  fit03aPd, # [df0] Agency ~ (Time | Name) + Time * Period
  fit04aPd, # [df0] Agency ~ (Time | Name) + Time * Period * Outcome
    show.reflvl = FALSE,
    show.intercept = TRUE,
    p.style = "numeric_stars",
    file = file.path(dir8, "summary-all-models-part-02-publication.html"))
tab0

	Agency			Agency			Agency			Agency
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	0.50 ***	0.49 – 0.50	<0.001	0.50 ***	0.49 – 0.50	<0.001	0.51 ***	0.50 – 0.51	<0.001	0.46 ***	0.44 – 0.47	<0.001
Time				-0.01 ***	-0.02 – -0.01	<0.001	-0.15 ***	-0.18 – -0.12	<0.001	-0.27 ***	-0.33 – -0.21	<0.001
Period [PB]							0.02 ***	0.02 – 0.03	<0.001	0.06 ***	0.04 – 0.07	<0.001
Period [PA]							0.02 ***	0.01 – 0.03	<0.001	-0.05 ***	-0.07 – -0.02	<0.001
Time × Period [PB]							0.19 ***	0.16 – 0.23	<0.001	0.33 ***	0.26 – 0.39	<0.001
Time × Period [PA]							0.08 ***	0.05 – 0.12	<0.001	0.25 ***	0.18 – 0.32	<0.001
Outcome [winner]										0.08 ***	0.06 – 0.10	<0.001
Time × Outcome [winner]										0.16 ***	0.08 – 0.23	<0.001
Period [PB] × Outcome [winner]										-0.05 ***	-0.07 – -0.04	<0.001
Period [PA] × Outcome [winner]										0.09 ***	0.06 – 0.11	<0.001
(Time × Period [PB]) × Outcome [winner]										-0.18 ***	-0.25 – -0.10	<0.001
(Time × Period [PA]) × Outcome [winner]										-0.21 ***	-0.29 – -0.13	<0.001
Random Effects
σ2	0.07			0.07			0.07			0.07
τ00	0.01 Name			0.01 Name			0.01 Name			0.00 Name
τ11				0.00 Name.Time			0.00 Name.Time			0.00 Name.Time
ρ01				0.18 Name			0.18 Name			-0.07 Name
ICC	0.08			0.10			0.10			0.08
N	870 Name			870 Name			870 Name			870 Name
Observations	169997			169997			169997			169997
Marginal R2 / Conditional R2	0.000 / 0.084			0.000 / 0.102			0.004 / 0.105			0.019 / 0.100
* p<0.05   ** p<0.01   *** p<0.001

# save.image(file=file.path(dir8, "session.RData"))
# load(file.path(dir8, "session.RData")

Report

result <- report::report(fit04aPd)

print(result)

We fitted a linear mixed model (estimated using REML and Nelder-Mead optimizer) to predict Agency with Time, Period and Outcome (formula: Agency ~ Time * Period * Outcome). The model included Time as random effects
(formula: ~Time | Name). The model's total explanatory power is weak (conditional R2 = 0.10) and the part related to the fixed effects alone (marginal R2) is of 0.02. The model's intercept, corresponding to Time = 0,
Period = PC and Outcome = loser, is at 0.46 (95% CI [0.44, 0.47], t(169981) = 61.71, p < .001). Within this model:

  - The effect of Time is statistically significant and negative (beta = -0.27, 95% CI [-0.33, -0.21], t(169981) = -8.35, p < .001; Std. beta = -0.57, 95% CI [-0.70, -0.44])
  - The effect of Period [PB] is statistically significant and positive (beta = 0.06, 95% CI [0.04, 0.07], t(169981) = 8.51, p < .001; Std. beta = 0.13, 95% CI [0.07, 0.20])
  - The effect of Period [PA] is statistically significant and negative (beta = -0.05, 95% CI [-0.07, -0.02], t(169981) = -4.05, p < .001; Std. beta = -0.23, 95% CI [-0.33, -0.14])
  - The effect of Outcome [winner] is statistically significant and positive (beta = 0.08, 95% CI [0.06, 0.10], t(169981) = 8.67, p < .001; Std. beta = 0.25, 95% CI [0.17, 0.33])
  - The effect of Time × Period [PB] is statistically significant and positive (beta = 0.33, 95% CI [0.26, 0.39], t(169981) = 10.07, p < .001; Std. beta = 0.69, 95% CI [0.56, 0.83])
  - The effect of Time × Period [PA] is statistically significant and positive (beta = 0.25, 95% CI [0.18, 0.32], t(169981) = 7.17, p < .001; Std. beta = 0.53, 95% CI [0.38, 0.67])
  - The effect of Time × Outcome [winner] is statistically significant and positive (beta = 0.16, 95% CI [0.08, 0.23], t(169981) = 4.16, p < .001; Std. beta = 0.33, 95% CI [0.18, 0.49])
  - The effect of Period [PB] × Outcome [winner] is statistically significant and negative (beta = -0.05, 95% CI [-0.07, -0.04], t(169981) = -6.58, p < .001; Std. beta = -0.15, 95% CI [-0.23, -0.08])
  - The effect of Period [PA] × Outcome [winner] is statistically significant and positive (beta = 0.09, 95% CI [0.06, 0.11], t(169981) = 6.65, p < .001; Std. beta = 0.37, 95% CI [0.26, 0.48])
  - The effect of (Time × Period [PB]) × Outcome [winner] is statistically significant and negative (beta = -0.18, 95% CI [-0.25, -0.10], t(169981) = -4.68, p < .001; Std. beta = -0.37, 95% CI [-0.53, -0.22])
  - The effect of (Time × Period [PA]) × Outcome [winner] is statistically significant and negative (beta = -0.21, 95% CI [-0.29, -0.13], t(169981) = -5.27, p < .001; Std. beta = -0.45, 95% CI [-0.62, -0.28])

Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation.